#version 430 core

#define LOC_X 16
#define LOC_Y 16

#define GRID_SIZE_X 64
#define GRID_SIZE_Y 64
#define UNIT_OFFSET_X 1.f/(GRID_SIZE_X+1)
#define UNIT_OFFSET_Y 1.f/(GRID_SIZE_Y+1)
#define GRAVITY vec2(0.0f, -98.0f)

layout(local_size_x=LOC_X, local_size_y=LOC_Y) in;

layout(binding=0) uniform sampler2D velocity;
layout(binding=1) uniform sampler2D levelSet;
layout(binding=0, rg32f) uniform writeonly image2D advectedVelocity;
layout(binding=1, r32f) uniform writeonly image2D advectedLevelSet;

uniform float timeStep;

void main(void)
{
	ivec2 gid = ivec2(gl_GlobalInvocationID.xy);
	if (gid.x <= GRID_SIZE_X && gid.y <= GRID_SIZE_Y)
	{
		// Semi Lagrangian lookup (Selle et al. 2007)
		vec2 vel_node = texelFetch(velocity, gid, 0).xy;
		vec2 pos_semilag = gid - timeStep*vel_node; // integer coords
		vec2 pos_semilag_tex = (pos_semilag+vec2(0.5,0.5)) * vec2(UNIT_OFFSET_X, UNIT_OFFSET_Y); // texture coordinates

		// Advect level set
		float level_semilag = texture(levelSet, pos_semilag_tex).x;
		imageStore(advectedLevelSet, gid, vec4(level_semilag, 0.f, 0.f, 1.f));

		// Modified MacCormack lookup for velocity
		vec2 vel_semilag = texture(velocity, pos_semilag_tex).xy;
		vec2 pos_mccorm_tex = pos_semilag_tex + vel_semilag * vec2(UNIT_OFFSET_X, UNIT_OFFSET_Y) * timeStep;
		vec2 vel_mccorm = texture(velocity, pos_mccorm_tex).xy;
		vec2 vel = vel_semilag + (vel_node - vel_mccorm)*0.5f;

		// Min/max of values contributing to Semi Langrangian scheme
		ivec2 pos_semilag_corner = ivec2(floor(pos_semilag));
		vec2 lower = texelFetch(velocity, pos_semilag_corner, 0).xy;
		vec2 upper = lower;

		vec2 lookup = texelFetch(velocity, pos_semilag_corner + ivec2(0,1), 0).xy;
		lower = min(lower, lookup);
		upper = max(upper, lookup);

		lookup = texelFetch(velocity, pos_semilag_corner + ivec2(1,0), 0).xy;
		lower = min(lower, lookup);
		upper = max(upper, lookup);

		lookup = texelFetch(velocity, pos_semilag_corner + ivec2(1,1), 0).xy;
		lower = min(lower, lookup);
		upper = max(upper, lookup);

		// Clamp advected velocity
		vel = clamp(vel, lower, upper);
		vel += GRAVITY*timeStep;
		imageStore(advectedVelocity, gid, vec4(vel, 0.f, 1.f));
	}
}